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NAVAL ARCHITECTURE 


Scale effects in cavitation experiments 
with marine propeller models 


Jan A. Szantyr 
Gdańsk University of Technology 


ABSTRACT 


The paper presents an overview of specific scale effects encountered in the cavitation experi- 
ments with marine propeller models. These scale effects result from unavoidable dissimilarity 
between model and full scale flow phenomena. They may influence first of all inception 
of different forms of cavitation, but they also are visible in development and desinence of 
cavitation on propeller blades. These scale effects may be divided into five main categories, 
which are described in detail. The influence of these five categories of scale effects on the 
different aspects of cavitation performance of marine propellers is discussed. 


Keywords: cavitation, marine propellers, model experiments, scale effect 


INTRODUCTION 


Model experiments with marine propellers are still an 
important element of the process of marine propeller design. 
A meaningful proportion of these experiments are devoted to 
cavitation phenomena themselves and to different consequences 
of cavitation. Marine propellers exhibit a variety of specific 
forms of cavitation, shown schematically in Fig.l. The main 
forms include sheet, bubble, vortex and cloud cavitation. All 
forms of cavitation are as a rule unsteady and during propeller 
blades motion they undergo a process of inception, develop- 
ment and desinence, either related to the unsteadiness of flow 
around the propeller or resulting from the inherent unsteadiness 
of the cavitation phenomena themselves. They may lead to 
many negative consequences, such as pulsating pressure field, 
hydro-acoustic signals, erosion etc. It is important that predic- 
tions of full scale cavitation behaviour of a propeller, made on 
the basis of model experiments, are realistic and accurate. This 
requires detailed knowledge of specific scale effects related to 
cavitation phenomena. 


Fig. 1. Different forms of cavitation encountered on marine propellers . 


With an exception of a few depressurized model basins, 
majority of model experiments with marine propellers are 
conducted in cavitation tunnels. As a rule these experiments 
are based on fulfilling two simple similarity criteria. The first 
criterion requires equality of model and full scale advance 
coefficients: 


V, V 
Jc cya (1) 
nsDs NyDy 
where : 
V — propeller advance velocity 
n — propeller number of revolutions 
D — propeller diameter 
index S — full scale (ship) values 


index M — model values. 


This criterion ensures similarity of the basic kinematics of 
model and full scale flow around the propeller blades. How- 
ever, this does not imply full similarity of flow in model and 
full scale, because Reynolds numbers in these two situations 
are markedly different. This leads to formation of different 
boundary layers in model and full scale. As a routine in the 
course of experiments it is only ensured that flow on the blades 
of a propeller model is supercritical, selecting sufficiently high 
number of propeller model revolutions to fulfil the following 
condition : 


_ 0.757n yD 


M 


Rey M&M >5.10° (2) 
where : 


c — blade section chord length 
at non-dimensional radius 0.75 
v — kinematic viscosity coefficient of water. 


Full scale propeller Reynolds numbers are in the range Re, = 
= 10’. In practice condition (2) does not always ensure fully 
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turbulent flow in model scale and regions of laminar flow may 
still be present on model propeller blades, especially near the 
hub. In full scale the flow is always turbulent from the leading 
edge. Other consequences of difference in model and full scale 
Reynolds numbers demonstrate themselves as differences in 
boundary layer separation or laminar/turbulent transition. These 
differences lead to the difference in open water characteristics, 
which require corrections [5] and they also have an influence on 
cavitation behaviour of the model and full scale propellers. 


The second criterion requires equality 
of model and full scale cavitation number : 


_ Pa tPs8hs—Pys _ o, = Prun=Pvm Q) 


> 2 Mi 2 
1 Ps VRS 4 Pu Vem 
where : 
Pp, — atmospheric pressure in full scale 
p — water density 
g — acceleration of gravity 
h — immersion of the point where cavitation number is 
defined 
py — critical pressure of water 
Prun ~ tunnel pressure at the point where cavitation number 
is defined 
V, — Teference velocity, either equal V or equal nD 


The second criterion is expected to ensure the similarity of 
inception, growth and desinence of different cavitation pheno- 
mena present on marine propeller blades. However, it is known 
for many years now that cavitation phenomena are far too 
complex to be scaled properly by a single simple parameter like 
the cavitation number. Fulfilling relation (3) leaves a number 
of important scale effects due to the following reasons : 


difference in static pressure distribution between the model 
and full scale propeller 

difference in water quality — first of all in nuclei content 
difference in boundary layer phenomena 

difference in time and velocity scales 

difference in size between model and full scale propeller. 


HHHMH 


These effects and their consequences are discussed in detail 
in the following sections. It is demonstrated that these scale 
effects may seriously affect the accuracy of full scale prediction 
of marine propeller cavitation performance made on the basis 
of model experiments. 


THE EFFECT OF THE STATIC 
PRESSURE DISTRIBUTION 


The typical scale coefficients of marine propeller models 
are in the range à = 20 — 30. Let us consider for example a ty- 
pical medium size propeller having the diameter D, = 6.0 [m] 
and its shaft centre submerged at 9.0 [m] below free surface. 
Its model may have the diameter D,, = 0.25 [m]. It is obvious 
that criterion (3) may be fulfilled exactly only at one selected 
depth. Usually either shaft centre or blade tip in upper posi- 
tion are selected for this purpose. This situation is explained 
in Fig.2, where an exemplary distributions of model and full 
scale cavitation numbers along propeller blades are calculated 
assuming 6, = 6,, at blade tip. 

As may be seen in Fig.2 there are quite meaningful differen- 
ces in the distribution of cavitation numbers along the blades 
in vertical position, resulting from the different static pressure 
distribution in model and full scale. These differences may 
result in discrepancies between model and full scale cavitation 
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6 =2.597 


FULL SCALE MODEL SCALE 


Fig. 2. Distribution of cavitation numbers along the blades 
of full scale and model propellers . 


pictures on the blades in vertical positions and also in the time 
history of unsteady cavitation associated with blade motion in 
the non-uniform ship wake velocity field. In order to provide 
an assessment of importance of this effect Fig.3 shows two 


Fig. 3. Comparison of cavitation extent observed at : 
a) o, = 6, (blade tip) and b) o, = c, (shaft centre). 


photographs of cavitation extent on a model propeller blade in 
two conditions, namely with criterion (3) fulfilled at blade tip 
and at shaft centre. As may be seen in this figure, the difference 
is quite substantial. The expected full scale cavitation extent 
at blade top position should be somewhere in between the two 
pictures shown in Fig.3. Of course as in most cases cavitation 
phenomena are concentrated near the blade tip, the picture in 
Fig.3a may be regarded as closer to full scale reality. However, 
if the sheet cavity extends along most of the blade leading edge 
or we are specifically interested in the hub vortex formation, 
these effects may be decisive in proper scaling of the results. 

The influence of static pressure distribution on the unsteady 
cavitation associated with blade rotation may be elucidated by 
comparing the variation of blade tip cavitation number with blade 
rotation in model and full scale, as shown in Fig.4. It may be seen 
in this figure that the rate of change of the cavitation number with 
blade position angle is much higher in full scale than in model 
scale. This may lead to a smaller angular extent of the cavitation 
phenomena on the blade in full scale than in model scale. 

o 


FULL SCALE 
3.200 3.200 
2.870 2.870 
2.597 
MODEL SCALE 
270 0 O[deg] 90 


Fig. 4. Variation of the tip cavitation number with blade position angle . 


0 = 315 [deg] 


o= 0 [deg] 


In order to demonstrate the importance ofthis effect, cavi- 
tation sketches for different blade positions taken at different, 
appropriately chosen cavitation numbers are shown in Fig.5. 
In this figure three blade positions are shown, namely 315 [deg], 
i.e. 45 [deg] before reaching top, 0 [deg] and 45 [deg] after pas- 
sing top. The upper row shows sketches at cavitation numbers 
corresponding to full scale at the respective blade positions, 
while the lower row shows sketches corresponding to model 
experiment with condition (3) fulfilled at blade tip in top position. 
These sketches should be of course treated only as a quasi-ste- 
ady approximation of the real processes, but still it gives some 
idea about possible differences in scaling to full scale. 

The scale effect discussed in this section may also influence 
other important cavitation characteristics of marine propellers, 
for example the propeller induced pressure pulsations, which 
are routinely measured in the cavitation tunnel and scaled for 
full scale. In this case setting condition (3) either at the blade 
tip or at the shaft centre leads to markedly different results 
which are shown in Tab.1. This table includes the amplitudes 
of the first three harmonics of pressure pulsations, given in 
[kPa] and starting with blade frequency. The pressure pulsations 
were measured in the cavitation tunnel and re-calculated into 
full scale values. The measurements refer to the point located 
immediately above propeller in the ship plane of symmetry. It 
may be seen that the alternative choice of setting condition (3) at 
blade tip or shaft centre produces almost 80 per cent difference, 
what correlates well with the pictures of cavitation in Fig.3. 
The true full scale prediction should be placed somewhere in 
between the two sets of values given in Tab.1, perhaps closer 
to the data for “blade tip” criterion. 


Tab. 1. Pressure pulsations amplitudes in[kPa] 


6,, = 6, blade tip | 6,,= 6, shaft centre 
1 
2 
3 


0 = 45 [deg] 


Fig. 5. Variation of sheet cavity with angular blade position observed in model scale (lower row) and simulated for full scale (upper row) . 
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THE EFFECT OF THE WATER QUALITY 
— NUCLEI CONTENT 


From the point of view of cavitation the water quality in the 
experimental installation or at sea is described primarily by the 
dissolved air content, nuclei content and its previous history. 
In general there are meaningful differences in these parameters 
between the laboratory and full scale situations, which are the 
source of specific scale effects. These effects affect first of all 
the inception of cavitation, but they also have some influence on 
further development and desinence of cavitation. The detailed 
mechanisms of these scale effects are described below. 

The classical cavitation similarity criterion (3) is based on 
the assumption that the critical pressure Pop, i.e. pressure at 
which the liquid starts to rupture under tension and cavitation 
develops, is equal to the vapour pressure of the liquid p, . How- 
ever, itis known for quite a long time that the amount of tension 
that a liquid may withstand is substantial and it depends first 
of all on the dissolved gas content, nuclei content and previous 
history of the liquid. All these three factors are strongly inter- 
-related with each other. Consequently, Pog in most practical 
situations is not equal p, and the cavitation number for inception 
of cavitation should be written in the following form : 


_ Prun=Pcr _ Prun=@v—Pr) _ 


Ol a 
TA ovè 


i (4) 
T 
xPVa 


=o+Ao 


where : 


tensile strength of the liquid 
inception cavitation number. 


Py 7 

Go, - 

Then the cavitation inception condition 
should take the following form : 


B C i = o (5) 
where : Cia is the non-dimensional minimum 
pressure coefficient on the tested object. 


It would be very convenient if the value of p, could be 
related directly to the nuclei content of the water. However, 
due to the complicated physics of cavitation inception this is 
not possible. Some guidance concerning this relation may be 
obtained from a simple condition for the stability of a single 
spherical vapour bubble, which is developed from the well- 
-known Rayleigh Plesset equation [1] : 


4S 
Pr=3R (6) 
where : 
S — surface tension of water 
R — bubble radius. 


Equation (6) shows that smaller bubbles can withstand 
higher tension and vice versa. In case of a random distribution 
of real nuclei the value of p, may be either negative or positi- 
ve, or in particular cases equal to zero. The liquid containing 
large bubbles exhibits negative tension strength in the sense 
that bubbles may grow due to diffusion without applying any 
tension to the liquid at all. On the contrary, liquid containing 
only a few small nuclei exhibits quite high positive tension 
strength. The dramatic effect ofp, on the cavitation phenomena 
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may be seen in Fig.6, showing cavitation phenomena on a wing 
at high positive p, , at p, = 0 and at negative p, . 


a) 


Fig. 6. Cavitation phenomena on a wing at high p, (a), 
at p,, = 0 (b) and at negative p, (c), quoted from [3] . 


As the value of p, depends strongly on the nuclei content of 
the water and also on the dissolved air content, these two para- 
meters should be controlled during experiments in the cavitation 
tunnel in such a way that consistent results are obtained and 
reliable scaling of the results into full scale is possible. Nuclei 
content in tunnel water may vary considerably, what is shown 
in Fig.7 taken from [1]. 

The influence of the nuclei distribution in water is visible 
first of all in determining the inception conditions for different 
forms of cavitation. This effect is illustrated in Fig.8, showing 
the variation in the cavitation inception number o, with nuclei 
concentration for three forms of cavitation: tip vortex cavita- 
tion, sheet cavitation and bubble cavitation. The effects shown 
in this figure may be explained in the following way : 


> bubble cavitation is generated by rather small tension di- 
stributed widely over the central part of the propeller blade. 
This tension is able to destabilize only relatively large and 
sparsely distributed nuclei, hence its strong dependence on 
nuclei concentration 


> sheet (surface) cavitation is generated by strong tension 
located in a narrow zone near the blade leading edge. This 
tension destabilizes even the smallest nuclei, hence its 
dependence on nuclei concentration is much weaker 


> tip vortex cavitation is initiated in a small volume around the 
centre of the vortex. Its inception depends on the probability 
of nuclei of suitable size flowing into that volume. This 
probability is the direct function of nuclei concentration, 
hence the dependence is strong. 


Untreated tap 
20 water 
Total with 
10 
T 
g9 : 
e Degassed tap 
z 20 water 
Q 
ka 
s 
a 
810} 
z 
[5] 
3 
d 


Filtered tap 


10 
Nuclei radius R [um] 


Fig. 7. Distribution of nuclei in water in three stages of degassing [1] . 


4 o[-] 


Propeller with 


bubble cavitation 

3 a 

2- 

1 4 

Tension 
OSS SSS 
Low Med. Med. High 
High conc. Low conc. 
4 o[] 

134 


Propeller with 
blade surface cavitation 


ne ei ëe O 


11- 
Tension 
10 T T T T 
Low Med. Med. High 
High conc. Low conc. 
toH 
20-4 
Propeller with 
45 tip vortex cavitation 
10-4 
Tension 
5 T T T T g 
Low Med. Med. High 
High conc. Low conc. 


Fig. 8. Influence of the nuclei concentration on the inception cavitation 
number for three forms of cavitation [2] . 


The influence of dissolved air on the cavitation inception 
condition may be qualitatively illustrated by Fig.9 showing 
experimentally obtained curves for inception of tip vortex, hub 
vortex, sheet (surface) and bubble cavitation for two concen- 
trations of dissolved air. 


poh 


Cavitation inception 


p 
C 
o 
Q 
9 
< 
= 
Pie 


Blade surface cavitation 


aaant 


Advance ratio, J [-] 


Fig. 9. Influence of dissolved air on cavitation inception 
on a marine propeller model [2] . 


The above described scale effect in cavitation experiments 
with marine propeller models is especially important for naval 
and other special applications. In these applications ensuring 
a completely cavitation-free operation is essential because of 
acoustic emission. It requires very precise knowledge about 
scale effects important for scaling cavitation inception condi- 
tion. This implies, among other things, information about nuclei 
content in sea water. Measurements of such nuclei content were 
conducted for example by the French Navy [2]. Results of the 
assessment of sea water tensile strength at different depth are 
shown for example in Fig.10. These results may serve as certain 
guidance in scaling of the results of model experiments with 
marine propellers into full scale. 


-100000 - 50000 0 


Pcr - Pv (N/m?) 


-150000 


Fig. 10. Results of measurements of sea water tensile strength 
at different depths [2] . 


In case of typical cargo ship propellers it is generally ac- 
cepted that routine guidelines for degassing water in cavitation 
tunnel and keeping the dissolved air content below about 30 per 
cent of saturation are sufficient to ensure the model cavitation 
test results corresponding to full scale situation. However, it 
is still important to know the nature of the scale effect due to 
the water quality. 
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THE EFFECT OF THE BOUNDARY 
LAYER PHENOMENA 


The boundary layer phenomena depend primarily on the 
value of the Reynolds number. This number is composed of the 
characteristic velocity, characteristic dimension and viscosity 
coefficient. However, in the case of cavitation scaling it may 
be distinguished between Reynolds number effect (strictly 
linked to the boundary layer phenomena) and separate effects 
of the flow velocity and object size, which have a completely 
different physical mechanism and are discussed separately in 
the following sections. 

In turn the “pure” Reynolds number effect may be divi- 
ded into (1) a difference in quasi-steady pressure distribution 
between model scale and full scale propeller blades and (2) 
unsteady velocity/pressure fluctuations inside the boundary 
layer. The first of these two effects is of lesser importance, 
but nonetheless it should be taken into account. The limited 
relative importance of this effect may be assessed on the basis 
of lift-related scale effect corrections for marine propellers as 
discussed for example in [7]. 

The second effect is much more important, particularly from 
the point of view of cavitation inception scaling. According to 
[6] pressure fluctuations inside the boundary layer on a blade 
section profile may reach quite high amplitudes, comparable 
with the mean value of pressure in the vicinity of the leading 
edge. These fluctuations may be added to the distribution of 
mean pressure along the blade section profile. 


ene 


MODEL SCALE Cc 


FULL SCALE a 


Fig. 11. Pressure fluctuations inside the boundary layer 
on a blade section in model scale and full scale . 


These fluctuations reach maximum amplitudes in the region 
of laminar/turbulent transition. Due to the difference between 
model and full scale Reynolds numbers this region of transition 
is located differently on the model scale blade section profile 
and in corresponding full scale situation (cf. Fig.11). In full 
scale this region is located close to the leading edge, where the 
point of minimum of the mean pressure coefficient C is also 
located if the profile operates at an angle of attack. In model 
scale this region is located further downstream, where the mean 
pressure on the suction side of the profile is already markedly 
higher. Consequently, at the value o, of the local cavitation 
number inception may take place in full scale, but not in model 
scale. Typical frequency of pressure/velocity oscillation in the 
transition zone is of order 1 [kHz], while the own oscillation 
frequency of cavitation nuclei between 20 and 50 microns (most 
frequently encountered in real flows) is of order 50 [kHz]. From 
this we may conclude that any of these micro-bubbles reaching 
the laminar/turbulent transition zone will be destabilized by the 
pressure pulsations and may cause inception of cavitation. If 
the local cavitation number is sufficiently low, say equal to o, , 
inception and development of cavitation will be experienced 
both in full and model scale. 
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THE EFFECT OF THE TIME 
AND VELOCITY SCALES 


It is generally accepted that inception of cavitation takes 
place when cavitation nuclei present in the flow are destabili- 
zed by the flow conditions and they undergo dynamic growth 
[5]. In order to achieve this, the nuclei must be acted upon by 
a sufficiently low pressure during a sufficiently long time. This 
situation is illustrated in Fig.12. In this figure C is the non- 

‘ è a n 4 P 
-dimensional pressure distribution along a blade section profile 
and o, means the theoretical cavitation inception number, based 
on the assumption Pop = Py - Due to the above described effect, 
the actual inception cavitation number will be lower and equal 
o It may be assumed that a certain value C of the multiple of 
mean pressure difference between local pressure and vapour 
pressure Ap and time At is required to initiate cavitation. This 
parameter C is inversely proportional to the typical size of 
cavitation nuclei (compare Eq. 6), so it may be assumed that 
its value is similar in model and full scale. Then the following 
relation may be developed: 


L X 
C=ApAt=ÊV2AC =A% = 
Pp 2 0 P Vo L 
X,/L (7) 
P Tal X 
=5 VoL fC- opal 
x,/L 
where : 
V, — characteristic velocity of flow at the blade section 
L — length of the blade section 
AC, — non-dimensional pressure difference (cf. Fig. 12). 
P 
oO 
V 
= 0; 
0 xX, X2 ”x 
L L L 


Fig. 12. The effect of time and pressure difference on cavitation inception . 


The following observations 
may be made from the form of equation (7) : 


> the same object tested at different velocities will have dif- 
ferent values of the actual inception cavitation number oj; 
this value will be smaller for smaller velocity; on the other 
hand it will be closer to o, for larger velocity 


> fortwo geometrically similar objects of different sizes tested 
at the same velocity the smaller object will exhibit smaller 
value ofo} 


> two geometrically similar objects of different sizes tested at 
the same value of V L will exhibit identical values of or. 


If the relation (7) is applied to the situation of the model 
and full scale marine propeller, it may be concluded that the 
characteristic velocity V, in full scale is about twice the model 
value and the characteristic length L in full scale is about 25 
times the model value. Consequently, the actual inception ca- 
vitation number ol will be visibly higher in full scale than in 
model scale. As a result, cavitation inception in full scale may 
happen earlier than indicated by the model experiments. This 
effect may be illustrated by the photographs shown in Fig.13. 
In these photographs cavitation on the so called Schiebe body 
is shown at two different flow velocities but at the same value 
of the cavitation number. 


Fig. 13. Schiebe body tested at the same cavitation number 
but at V, = 8 [m/s] (a) and at V, = 14 [m/s] (b), quoted from [3] . 


THE EFFECT OF THE OBJECT SIZE 


One could assume that if the scale effect due to the differen- 
ce between model and full scale Reynolds number is properly 
taken into account, then there is no need to analyse separately 
the effect of object size. However, the complicated nature of 
cavitation leads to a separate scale effect of the object size. This 
may be justified by the following reasoning : 


= two geometrically similar objects of different sizes will 
exhibit similar regions of low pressure having size propor- 
tional to the object size 


= ifthe same distribution of the cavitation nuclei in the onco- 
ming flow is assumed in both cases, then more nuclei will 
enter the low pressure region around the larger object 


= consequently, the probability of destabilization of nuclei 
and cavitation inception is higher for a larger object. 


This effect may be illustrated by the photographs shown in 
Fig.14. However, it should be kept in mind that in the laboratory 
experiments there are always many different factors influencing 
the results at the same time. The results shown in Fig.14, but 
also in Fig.13 may be the joint consequence of effects mentio- 
ned in this section and in two preceding sections together. 


Fig. 14. Three geometrically similar Schiebe bodies 
of different sizes tested at the same cavitation number : 
a) diameter 15 [mm], b) diameter 30 [mm], 

c) diameter 60 [mm], quoted from [3]. 


CONCLUSION 


The above presented analysis of the specific scale effects, 
influencing full scale predictions of the marine propellers 
cavitation behaviour made on the basis of model experiments, 
leads to the following conclusions : 


O most of the specific scale effects are more important for 
inception of cavitation on marine propellers and less 
important for the fully developed cavitation phenomena; 
consequently, these scale effects are more important for 
design of “silent propellers” (naval, research, passenger 
ships etc.) than for routine cargo ships propellers 


O the unavoidable effect of static pressure distribution may 
affect both inception and development of cavitation; it 
may also generate meaningful differences in the propeller 
induced pressure pulses; the influence of this scale effect 
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may be minimized by proper selection of the cavitation 
number in model experiments 

O the water quality effect, in particular the distribution of ca- 
vitation nuclei in the experimental facility and at sea, is the 
source of important scale effect on cavitation inception; at 
present it can not be fully taken into account during scaling 
of the model results due to the lack of information about 
true full scale situation 

O the boundary layer scale effect may lead to meaningful 
under-prediction of the cavitation inception condition in 
full scale 

O the velocity and time scale effect may also lead to under- 
-prediction of the full scale cavitation inception condition 

O the effect of object size also points to a possibility of under- 
-prediction of the full scale cavitation inception condition 

O consequently, the results of model experiments performed 
for prediction of cavitation inception on full scale propellers 
must be taken with special care and the designers should be 
aware of the fact that all specific scale effects cause the full 
scale propeller cavitation behaviour visibly less optimistic 
than indicated by model experiments. 

NOMENCLATURE 

c - blade section chord length 

C -non-dimensional pressure coefficient 

C in - minimum value of the non-dimensional pressure coefficient 

AC, -non-dimensional pressure difference 

D -propeller diameter 

g - acceleration of gravity 

h - immersion of point where cavitation number is defined 

J - propeller advance coefficient 

L - characteristic length 

n - propeller number of revolutions 

p, -atmospheric pressure 

Ap -difference between local pressure and vapour pressure 

Per ~ Critical pressure 

p, -cavitation inception pressure 

p, -tensile strength of the liquid 

Prun ~ Pressure inside cavitation tunnel 

Py -vapour pressure 

R -cavitation bubble radius 

Re = - Reynolds number 

S - surface tension of water 

At -time necessary for cavitation inception 

V -propeller advance velocity 

Və -characteristic velocity 

V, -reference velocity 

> - non - dimensional coordinate along profile chord 

0 - propeller blade position angle 

A - model scale coefficient 

v - kinematic viscosity coefficient 

p - density of water 

(o - cavitation number 

Ao -change in cavitation number due to tensile strenght of water 

6,  - inception cavitation number 

ol - inception cavitation number including water tensile strength 
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Miscellanea 


Every year the Polish Society of Naval Architects and 
Marine Engineers, TOP KORAB, traditionally arranges, 
apart from its current activity, open meetings during which 
are discussed topics in advance announced dealing with 
shipbuilding and maritime economy. 


In the year 2005 the following problems were dealt with 
by Gdansk Division of the Society : 


* Prospects of Polish shipbuilding industry in taking 
advantage from the current boom in the worldwide 
market 

Drilling unit equipment made in Poland 

Fuel system blocks for ship diesel engines 

The fastest-in-the-world computers and their techno- 
logical trends 

To day and to morrow of the firm Hydrobudowa 
S.A. 


And in Szezecin Division of the Society : 


+ ot 


Optimization models for engineering and economy 
Organization of the Conference on : Past, present and 
future state of Polish shipbuilding industry 


++ 


In 2005 the most important event for the Society was 
the signing of the agreement on establishment of RINA/ 
KORAB Joint Branch. The agreement was signed on 
behalf of TOP KORAB by B. Banach, its Chairman and 
J. Czuczman, its Vice-chairman, and on behalf of RINA 
by T. Blakeley, RINA Chairman. The document was also 
signed by Prof. J. W. Doerffer, Honorary Member of TOP 
KORAB, Fellow Member of RINA. 
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Effective fatigue damage summation 


Janusz Kolenda 
Gdansk University of Technology 
Polish Naval University 


ABSTRACT 


The cumulative fatigue damage caused by uniaxial variable stress is considered. To facilitate 
the lifetime prediction at complex stress patterns, has been defined an effective stress of 
Y constant parameters, claimed to be equivalent to the sequence of original stresses in terms 
of fatigue life of material. The equivalence conditions are based on the theory of energy 
transformation systems. To the effective stress the Palmgren-Miner rule is applied. Thereby 

the cycle counting is not required. 


Keywords : fatigue, design criteria, damage cumulation, stress modeling 


INTRODUCTION 


Most engineering components are subjected to variable 
loading. The fatigue damage in materials under cyclic loads 
of variable level, direction and/or mode, is accumulated at an 
unsteady rate. Though many cumulative damage theories have 
been developed, none of them enjoys universal acceptance. 
Each fatigue damage model can only account for one or several 
phenomenological factors, such as load dependence, multiple 
damage stages, non-linear damage evolution, load sequence 
and interaction effects, mean stress, overload effects and small 
amplitude cycles below fatigue limit. The most convenient is 
the Palmgren-Miner rule [1,2], and it has been adopted and 
used extensively, mainly because of its simplicity. According 
to the original concept : 


v the fatigue process is cumulative, 

v the fatigue effect is proportional to the work done by active 
loads, 

v the increment of damage caused by n stress cycles of 
constant amplitude can be estimated as n/N, where N is 
the number of stress cycles which would cause failure of 
material in the same load state. 


As to the loads of variable amplitude, the concept holds that 
the total damage can be estimated as the sum of the damage 
increments, each corresponding to a specified stress level. This 
can be symbolically expressed by the equation : 


pay (1) 
a Nj 
where : 


D - total fatigue damage which equals unity at the failure 
n, - number of stress cycles at j-th stress level 

j- number of cycles to failure at j-th stress level 
w - number of stress levels. 


The quantity N, can be calculated as follows [3] : 


N= ga -Za = 


2j 
where : 


Ww 


- fatigue strength coefficient 
(in equation of relevant S-N curve) 
- fatigue strength exponent 
- yield strength - for ductile materials, 
ultimate strength - for brittle materials 
0, - stress amplitude at j-th stress level 
o , - mean stress at j-th stress level. 


mj 


As 


Similar approach is applied to loads of variable level, 
direction and mode, namely : 


wh nj f 
DEDA. >; 1=X, y, Z, Xy, yZ, ZX (3) 
j=l + \ij 


1 


where the indices i = x,y,z are associated with normal stresses, 
and i = xy, yz, zx with shear stresses in Cartesian coordinate 
system. 

The Palmgren-Miner rule implicitly assumes that the damage 
cumulation is independent of the order in which the stress cycles 
of different levels are applied and of the order of load modes. 
The assumptions have been found inadequate in certain cases 
[4]. Moreover, the cycle counting may be a source of uncer- 
tainties and D values may differ significantly from unity at the 
failure. Therefore several other cumulative damage hypotheses 
have been proposed [5] but their applicability varies from case 
to case. In this context a relatively simple approach based on 
the Palmgren-Miner rule and effective stress models, can be 
mentioned. Such approach for zero mean variable stresses has 
been presented in [6]; it will be used below for nonzero mean 
variable stresses. 
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EFFECTIVE STRESS 


In accordance with the theory of energy transformation 
systems [7] the following can be stated. 

An effective stress model and an original stress can be re- 
garded equivalent in terms of fatigue life of material if during 
its service life the internally and externally dissipated energies 
under effective and original stresses are equal, respectively. 

To illustrate the determination of the effective stress, let 
us consider a variable stress of i-th mode in the high-cycle 
fatigue regime : 


G,(t)=6,,+6,(t) ; j=1,2,..., W (4) 
where c(t) is a zero mean time-varying stress. 
The sequence of stresses (4) can be modeled 
by the effective stress : 
S(t) = One + Oe SIN®,t (5) 


where ©» 6, and œ, are its mean value, amplitude and circular 
frequency, respectively. The equivalence conditions correspon- 
ding to the above a statement are [6] : 
t; 
J 


faoa- ava (6) 
j=l 
tj ; 
fa pde > [aja (7) 
0 j=l 
where : 


t- the duration of j-th load state, dot denotes differentiation 
with respect to time, and T : 


Tye 
j=l 


Eqs (6) and (7) determine parameters of the effective stress. 
For example, when : 


(8) 


S(t) =0,,,+ 0; Sin@ jt TERZ aw (9) 
these equations yield : 
1l w j 1/2 
Ome= TAPA 
j=l 
" 1/2 
O.= FoI] (10) 
jel 
x 22 
2 tjojoj e 
H 
O. = = 


In Eq (9) 5; is the stress amplitude 
and ©- the stress circular frequency at j-th level. 


FATIGUE DAMAGE 


On determination of the effective stress parameters and the 
number of stress cycles during w load states of i-th mode : 


e T (11) 
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the following is postulated : 


> the fatigue process at the effective stress is cumulative, 

> the fatigue effect is proportional to the dissipative energy, 

> the increment of damage, corresponding to n, cycles of the 
effective stress, can be estimated as follows : 


D.= N. =D, (12) 
where 
-K [,_ Gme\" 
N.= or 1 R | (13) 


Additionally, if the load is variable then the total fatigue 
damage based on the concept of the effective stress 
is assumed to be : 


D, =ġ Dai 


Of course, with iie presented approach the aforementioned 
drawback of the Palmgren-Miner rule, associated with ignoring 
the load sequence effects, cannot be avoided, but the cycle 
counting is not required. 


EXAMPLE 
Task : 


“* Compare the cumulative fatigue damages D and D. ifthe load is 
of i-th mode and the stress sequence consists of n, cycles 
of the stress : 


=E Yes (14) 


6,(t)=6,,+6,sin@,t 
and n, cycles of the stress : 
S(t) = Ow + 6, Sin Wt 

in the high-cycle fatigue regime when : 


mj j _ 
RPF?! pty : i=l? 
where : 


F - the fatigue limit under fully reversed stress of i-th 
mode 

L - the maximum stress amplitude o satisfying equation 
of i-th S-N curve 


No™=K 
above which the low-cycle fatigue may occur. 
** Consider the following cases : 
Ny=N,=N ; 6,,/R=0.5 ; o,,/R=0.5y 
0, =0,=90 ; ®, =, =0 


MAMAN 5 Oy = 


B m. m 
02= 70; ; Oy =0,=0 

n=i=n0 ; Ogi Ope Om 
> ye WE 


Solution : 


A. According to the Palmgren-Miner rule and Eq.(2), 
the fatigue damage D equals : 

n n 

a r ay ee Z] 


PEN, N, N, N, 


where : 
Ni=—rl-0.5)" No =n l- 0.57)" 


The parameters of the effective stress are calculated from 
Eqs (10) for t, = t,= 0.5T, as follows : 


Oo 1 o? o? ne 5 1/2 
=| k mEt, ‘al = [0.1250 +ô) 


TUR R? 
1/2 
=O 


1 2 2 
0. = m +t,05) 


2 241/2 


22 
[so +1,0505 
O.= 


2 2 
t,o, +t,0, 


Since the load duration is : 


T ai 4 20. n = 2@n 
©] 0o “ @ 
the number of cycles of the effective stress amounts to : 
Oe 
n,.==—T=2n 
e 2m 


and the number of its cycles to failure is : 


N.= 5° [1- Sue) K fi -[o.125 aay] | 
Hence 
D.= ene 1/2)m 
K{-[o.12sa+y?] | 
and 


D i i 
sisl- 4, 
Dz lage *@-05y) | 


a (15) 
{10.125 0+y7)] "| 


The results of the exemplary calculation of the ratio D/D, 
are given in Tab. 1. 


Tab. 1. Results of calculation with the use of Eq. (15) . 


The parameters of the effective stress are : 
1/2 
One=[4(0.5TS3, +0.5TS3,)| = 6, 
1 1/2 1/2 
o= kosto? +0.5T7?0) =|0.5( n ’)| o; 
7 | 0.5T@o? +0.5Tw?y207 |” 
i 0.5To? + 0.5Ty’o7 
so that : 
K Gal 
= 7 1- 
[0.50 +y] or R 
and 
MEA 
n No Nwo m 
D _ N, N, = l+y (16) 
D m 2n 7 2\] m/2 
: a 2\0.504 y?)| 


e 


Exemplary values of D/D, are presented in Tab. 2. 
Tab. 2. Results of calculation with the use of Eq. (16) . 


C. In Tab. 3 the values of the ratio D/D, are presented as 
calculated for the stress of variable frequency at different 
relative frequency ranges : 


A®@ _ | - y| 
On VY 
Here 


Aw=|o,-@,| 


is the difference between higher and lower circular 
frequency of the stresses ©; (t) and o,(t), and: 


Om= y 00 


is the central circular frequency of the stress cycles. 


For the period : 

2m 2m 2mn 1 
Ta tytty= im +n a (+g 

1 2 


the parameters of the effective stress are : 


Ome= Om > 0,=90 > O.= yo 
and the number of its cycles is : 


By applying Eqs (2) and (13), the number of cycles to failure 
under considered stresses is the same, namely : 
K 


O m 
N= N,= N.= sr (I-52 


m 


POLISH MARITIME RESEARCH, No 4/2006 13 


MARINE ENGINEERING 


MARINE ENGINEERING 


Consequently : 
D _2n_ A a7) 
D. n, l+y 


Tab. 3. Results of calculation with the use of Eq. (17) . 


From Tab. 1, 2 and 3 it can be observed that the linear sum- 
mation of the fatigue damage corresponding to the effective 
stress, based on the theory of energy transformation systems, 
leads to less conservative results for the stresses of variable 
level (in particular at higher values of the fatigue strength expo- 
nent), and to even more conservative results for the stresses of 
variable frequency (especially at the relative frequency ranges 
greater than unity). The observed discrepancies between the 
values of the fatigue damages D and D, result from different 
assumptions concerning the work of the active loads in the 
Palmgren-Miner rule and the dissipative energy in the equiva- 
lence conditions (6) and (7). 


CONCLUSIONS 


O For majority of engineering components, fatigue process 
under service conditions involves variable loading history. 
The cycle counting technique and the Palmgren-Miner rule 
are then frequently used for fatigue lifetime prediction. 
However, at complex stress patterns the cycle counting 
procedure may be complicated. 


O The effective stress defined in the foregoing text, makes 
it possible to avoid this difficulty if the stress process is 
integrable in respect of time over the service period. 


Q From the exemplary calculations it follows that the li- 
near summation of fatigue damages corresponding to the 
effective stress leads to less conservative results than that 
corresponding to the original stress of variable level , and 
to even more conservative results when the frequency of 
the original stress is variable. 


NOMENCLATURE 

D - fatigue damage 

D, - fatigue damage corresponding to the effective stress 

D„ - fatigue damage corresponding to the effective stress of i-th 


mode (i=x,y,Z, Xy, YZ, ZX) 

- fatigue limit under fully reversed stress 

fatigue strength coefficient 

- maximum stress amplitude satisfying equation of the S-N 
curve above which the low-cycle fatigue may occur 

- fatigue strength exponent 

- number of stress cycles 

- number of the effective stress cycles 

number of stress cycles at j-th level 

- number of stress cycles of i-th mode at j-th level 

- number of zero mean stress cycles 

- number of the effective stress cycles to failure 

- number of stress cycles to failure at j-th level 


ZPUP PPB 


AZ 


— 
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number of stress cycles to failure at i-th mode and j-th level 

- yield strength - for ductile materials, ultimate strength - for 
brittle materials 

- time 

- duration of j-th load state 

- duration of w load states 

- number of stress levels 

number of levels of i-th mode stress 

- multiplying factor 

- stress amplitude 

- amplitude of the effective stress 


att 


aasan 


o 


effective stress 7 
o; - amplitude of the stress ©; (t) in Eq. (9) 


o((t) - zero mean time-varying stress of j-th level 

9; (t)- variable stress of j-th level 

6,, - mean value of the stress 

Oe -~ mean value of the effective stress 

0, -~ mean value of the stress 5; (t) 

œ  - circular frequency 

®, - circular frequency of the effective stress 

œ, - circular frequency of the stress 5; (t) in Eq. (9) 
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Processing methods of measured runs 
of torque and angular speed 
of combustion engine propulsion system 
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ABSTRACT 


In this paper four processing methods of cyclic data runs, namely : Fourier interpolation 
and transformation, least-squares approximation by means of trigonometric polynomials , 
synchronous averaging and movable approximation , were presented by using, as an example, 
discretely measured runs of torque and angular speed. It was demonstrated that for extrac- 
ting a useful signal out of disturbance background and decomposing the disturbances the 
appropriate methods are the synchronous averaging and multiple movable approximation. 
To emerge significant spectral lines from a spectrum the least-squares approximation by 


using trigonometric polynomials can be applied. The method contains a criterion for signal filtration and 

is insensitive to disturbances, run truncation and sampling irregularity. Fast Fourier Transform (FFT) may 

have an auxiliary significance for preliminary spectral analysis. Its main disadvantage is that in order to 

build a smoothing-out filter its parameters should be determined in advance, e.g. by applying one of the 
remaining methods. 


Keywords : Ship piston engine, torque, angular velocity, data processing, 
Fast Fourier Transform (FFT), least-squares approximation by using trigonometric polynomials (TLSA), 
movable approximation, synchronous averaging. 


INTRODUCTION 


Aconvenient source of diagnostic information about work 
of ship combustion engine cylinders may be the torque measu- 
red on shafts of propulsion systems and angular speed measured 
by using proximity detectors placed over flywheels and free 
ends of shafts or other elements kinematically connected with 
engine crankshaft. In runs of the above mentioned quantities an 
information about work of all cylinders and entire propulsion 
system is aggregated that is why they are interesting from the 
point of view of possible gaining a diagnostic information from 
them. On the torque and angular speed runs natural vibrations 
of propulsion system and its elements, variations of consumer’s 
power demand and technical state changes of propulsion system 
are superimposed. Runs of the above specified quantities are 
characterized by periodicity which results from cyclic work 
of piston engines. Sampling is timed by pulses coming from 
angular position of crankshaft. 

Searching for different symptoms in measurement series 
including those of technical state changes requires applying ma- 
thematical tools for making it possible to extract useful signals 
out of the background of other signals and disturbances. 


INTERPOLATION AND APPROXIMATION 
OF RUNS BY MEANS 
OF TRIGONOMETRIC POLYNOMIALS 


Essence of Fourier interpolation 
and transformation 


In order to highlight differences in results of data processing 
realized by using other methods it seems reasonable to remind 
in this work that the crucial element of Fourier interpolation is 
expansion of the periodic function with the period T : 


f(t)=f(tt+kT) k=1,2,3,.. 
into the infinite trigonometric series of the form : 
a = ; 
f(t) = T + > (ax cos 2mVv,T+b, sin2nv,t) (1) 
k=l 
where : v, = k/T — frequency. 


The constant a, and coefficients a, and b, are determined by 
means of commonly known integral formulae. Dimension of 
the argument t results from physical phenomena represented 
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by a measured run and is usually considered identical with 
time or axis of natural numbers ordering the set of samples of 
measurement series. 


Eq. (1) can be expressed in the following form : 


oo 
f(t)= Ao + >, A,cos(2nv,t- ®,) (2) 

k=l 

a b 
Aye ; A, =a? + b? ; O,.=arctg a (3) 
k 
where : 
A, — constant (mean value of series) 
A,» ®, — amplitude and phase of n-th spectral line, respec- 
tively, [2]. 


Ordered runs of the quantities A, and ® constitute spectra 
(images) of amplitudes and phases of a considered run in fre- 
quency domain and are Fourier transformation of the analyzed 
run. Initial run can be reconstructed by means of Fourier inverse 
transformation. Filtration in frequency domain is realized by 
means of relevant operations on spectra e.g. by multiplying 
a given spectrum by suitably selected weighting functions prior 
to making the inverse transformation. In particular if some 
spectral lines or their groups are removed it will constitute 
a frequency filtration. 

Generally the function f(t) and its image are considered in 
complex domain. Transition to the complex form is based on 
the known Euler formula [2]. In this work to apply the complex 
description was not necessary. However its knowledge is neces- 
sary in the case of making use of a software for determining 
transforms where results are given only in the complex form, 
e.g. as in Excel software. 

The Discrete Fourier Transform (DFT) concerns the 
transformation of discrete data series described with the finite 
constant sampling period At where form of the sampled func- 
tion is no longer known, that means that to shorten sampling 
period, if required, is not possible. If a series is finite (i.e. 
it contains N samples) then the series corresponds with the 
measurement series. 

The idea of transformation maintains unchanged and in this 
case it consists in interpolation of series by means of trigono- 
metric polynomials, as shown above. 

In connection with the transition from continuous function 
to its discrete representation, two problems arose already in 
theoretical phase: of Nyquist frequency and of aliasing [1, 2]. 

In 1918 H. Nyquist announced that 2v,T samples are suffi- 
cient for approximate representing the signal whose spectrum is 
limited by the frequency v, , and duration time is T. In 1949 K. C. 
E. Shannon published the proof of the theorem on run sampling 
and that on run reconstruction [2]. The Shannon theorem on run 
sampling cannot be expressed in a simple way and it is better to 
make use of the Shannon theorem on run reconstruction, which 
states that if for the sampling frequency v the relationship v > Vv, 
is fulfilled then the function f(t) is unambiguously reproducible 
on the basis of its discrete values [2]. In practice are applied 
greater sampling frequencies depending on sampling purpose 
[2] since for reproduction of a given function a sufficiently 
complex Shannon’s interpolation polynomial should be used 
and that the sampling process parameters do not satisfy the 
assumptions of Dirac’s pulse comb function. 

When a run whose spectrum contains spectral lines of 
frequencies greater than v, is sampled with the frequency v, 
then the aliasing phenomenon will occur [1, 2]. In order to 
be protected against this phenomenon, during measurements 
low-pass filters or, if possible, a sampling frequency increased 
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up to a safe value, is applied. As possible recording is limited 
only to finite representations of the investigated measurement 
runs the next problem arises, namely that of leak which consists 
in appearing, in a given spectrum, not existing spectral lines 
resulting from cutting-off sinusoidal components. Therefore it 
follows that the interval of measured series does not constitute 
an integral multiple of every sinusoidal component. An entire 
branch of mathematics and measurement engineering has 
been developed to elaborate suitable time filters (windows) to 
mitigate influence of the phenomenon on derived frequency 
spectrum [1, 2]. 

In the domain of practical computer-based spectral analyses 
the basic notion is Fast Fourier Transform (FFT). From the point 
of view of transformation theory the essence the FFT does not 
differ from that of the DFT. All the above given comments 
dealing with the DFT remain valid. In some elaborations any 
distinction between the DFT and FFT is not introduced [2]. The 
FFT problem consists in the computation algorithms owing 
to which it has been possible to greatly decrease number of 
main computational operations. Many kinds of FFT software 
have been developed. Cooley and Tukey, who presented their 
algorithm [6] in 1965, are regarded as the authors of the com- 
puter-based FFT. 

The FFT software is offered in many recognized program 
packets such as e.g. Matlab, Mathematica, Excel. They may 
differ from each other in different algorithms for determining 
the FFT , and results obtained by means of different programs 
may appear not identical for some parameters of the same 
data sets. 


Analysis of results of measurements of torque 
and angular speed runs with the use of the FFT 


Transforms of example measurement series and several 
its testing modifications were determined in order to compare 
features of the analysis of measured runs made by using the FFT 
and least-squares approximation with application of trigonome- 
tric polynomials. To this end was used the torque and angular 
speed runs (shown in Fig.1) measured on a real object. 


ia WI 
CAM 
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Fig. 1. Measured runs obtained by means of a photo-electric 
speed indicator : M — torque (dimensionless), 
op [rad/s x 20] — angular speed pulses related to its mean value . 


The above presented measurement results were obtained by 
means of the experimental unit of photo-electric disc torque 
meter of ENAMOR [4], on the test bed of one-cylinder two- 
-stroke combustion engine [3]. From every shaft rotation 30 
samples of both torque and angular speed were derived in the 
domain of shaft rotation angle, that results from the torque 
meter construction. Values of the torque M were left dimen- 
sionless as it does not matter for discussion on run processing. 
This example is interesting since characteristic symptoms and 
disturbances are distinctly present in the measured runs, espe- 
cially those of the torque M. 

The amplitude spectrum of the run of the torque M, de- 
termined by means of the FFT nested in Matlab, is shown in 
Fig.2. 
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Fig. 2. The amplitude spectrum (FM) of the run of the torque M (Fig.1) 
derived by using Matlab software; M, — amplitude of spectral lines, 
v— frequency, k — FFT natural axis (numbers of spectral lines) . 


Majority of FFT programs (among them also Matlab and 
Excel) take measurement input data as complex numbers. In 
such case the number of symmetrical spectral lines N = 2K 
where N — number of samples of measurement series, and a con- 
stant is calculated. In the considered case where N = 480, apart 
from the constant , just 480 values of amplitudes and phases, 
240 out of which constitute real values, are calculated. There are 
programs where only real values are determined. Phase spectra 
are not presented in this work as they are deemed useless for 
presentation and discussion of the problems in question. 

To illustrate the symptoms which appear in results of real 
measurement data processing with the use of the FFT, was 
determined the test run M, whose amplitude spectrum QM, 
is shown in Fig.3. 
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Fig. 3. The amplitude spectrum QM, of the test run M,, determined 
by applying the least-squares approximation method to the run M (Tab.1) : 
M,=3.1 ; M,= 10.46. 


The spectral lines of the spectrum (Fig.1) were obtained 
by using the least-squares approximation of the run M (Tab.1). 
In the figures spectral lines for v > 9 as well as those of small 
amplitude values were omitted in order to obtain better quality 
of imaging the spectra. Only 10 spectral lines were left due to 
which it was possible to better demonstrate influence of some 
factors on transformation results. 


The amplitude FFT of the M, test run, 
achieved by using the Matlab software is shown in Fig.4. 
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Fig. 4. The amplitude spectrum FM, (FFT) of the M test run, 
obtained by means of the Matlab software . 


As compared with the spectrum QM, (Fig.3) in the spectrum 
FM, (Fig.4) characteristic side spectral lines (stripes) appeared 
at non-integral frequencies. 

The next test run, M. , was formed by removing 10 samples 
from each its end. The M. run contains N = 460 samples. 

The FFT (obtained by using Matlab) is characterized by 
a much greater number of side spectral lines (stripes) especially 
in the range of low frequencies (Fig.5). 
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Fig. 5. The amplitude spectrum FM, of the M_, test run, obtained 
by means of Matlab software; M,, — the run formed from M,,run 
by removing 10 samples from each its end . 


The phenomenon may constitute an obstacle in searching 
for characteristic spectral lines of lower amplitudes. Such cases 
can be met in measuring practice. Runs recorded without syn- 
chronizing and controlling pulses of starting-up and stopping 
the measurements will be generally truncated and even in the 
case of synchronization the truncation effect may happen for 
some sinusoids. 

In the case of application of some programs, removal of 
a part of samples or extension of measurement series by adding 
samples of zero values may be necessary. For instance the FFT 
program nested in Excel software demands sets of samples of 
size N = 2', where | - natural number. 


In Fig.6 is shown the transform FM., of the test set M_., 
obtained from the set M, of N = 480 samples by complementing 
it with zeros up to N = 512 samples. 
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Fig. 6. The amplitude spectrum FM of the test run M, formed 
by applying Excel software: M,, — test run obtained from the set M, 
of N = 480 samples by complementing it with zeros up to N = 512 samples . 


In the spectrum FM., (Fig.6) the side stripes of significant 
values can be observed. The disturbances are analogous to those 
as in the case of truncation of the run (Fig.5). 


The successive test run, M_,, was obtained by adding , to the 
run M, , pseudo-random disturbances of sampling frequencies 
generated by computer from the interval <-5 , 5> (of uniform 
distribution). As observed in Fig.7 the spectrum FM, contains 
a significant number of spectral disturbing lines. 
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Fig. 7. The amplitude spectrum FM „ of the test run M,, formed 
by applying Matlab software: M_, — test run obtained from the run M, 
by adding pseudo-random disturbances from the interval <-5 , 5> . 


In the example shown in Fig.7 the spectral lines resulting 
from disturbances can be considered significant. 
To prevent occurring such symptoms during measurements, 
as mentioned above, low-pass filters, synchronized recording 
and various time windows are applied [1, 2] . 


POLISH MARITIME RESEARCH, No 4/2006 17 


MARINE ENGINEERING 


MARINE ENGINEERING 


Approximation of runs by means 
of trigonometric polynomials 
and least- squares method (TLSA) 


It should be stressed once more that any transform obtained 
by using the FFT contains complete measured information, i.e. 
together with disturbances of analyzed time series, hence this 
is time series transformation into domain of frequency and 
phase, performed in the way of interpolation of measurement 
data by means of trigonometric polynomials of a given uniform 
(natural) frequency scale. 

Determination of spectral lines having their parameters 
close to real ones, which is connected with elimination of di- 
sturbances, can be realized by using the least-squares method, 
in the simplest case by assuming a power polynomial, as the 
orthogonal base, and sum of pairs of terms of a trigonometric 
polynomial : 


U K 
y=} ct" + > (a, sin v,t +b, cos vt) 
ial 


u=0 
p (4) 
t=2n— ; U+2K<N 
T 
where : 
C, — power polynomial coefficients 
a,,b, — trigonometric polynomial coefficients 
k — spectral line number 


t — angular axis, rad 

v — dimensionless spectral line frequency, v, > 0 

T — measurement axis (time, rotation angle, measure- 
ment point’s number) 

— period of main measurement axis 

— number of power polynomial coefficients 

number of trigonometric polynomial co- 

efficients 

— size of set of measurement samples. 
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The power polynomial has been introduced to approximate 
the constant and trends. In the Matlab software linear trend of 
data can be removed by using the command ,,detrend”. The 
trend is subtracted off the run by means of least-squares linear 
approximation. In the model (4) the power polynomial can be 
replaced by other functions t, and varying frequency can be 
introduced, that may be reasonable in the case of analyzing 
some runs of specific models. However, level of computation 
outlays may be limiting in such case. 

The coefficients c, and a, , b, are determined from the 
condition for minimum of sum of squares of deviations, and 
this sum is used to assess approximation quality and extract 
significant spectral lines. Values of amplitudes and phases are 
determined by using the formulae (3). The model (4) differs 
from the interpolation formula (2) by that values of frequen- 
cies are not to be multiples of the fundamental period and total 
number of determined coefficients and spectral lines must be 
smaller than that of samples, N, to ensure some number of 
degrees of freedom necessary to make it possible to perform 
statistical significance assessment of determined coefficients 
and spectral lines. 

The significant spectral lines are included one-by-one or 
collectively into the model with taking into account a degree 
of decreasing the sum of squares of deviations after their inclu- 
sion. If the successively included spectral lines do not make 
the sum of squares of deviations to decrease significantly then 
the approximation process will be stopped. 

In this work to analyze the torque and angular speed 
runs (Fig.1) the least-squares approximation algorithm using 
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trigonometric polynomials (the TLSA method) in accordance 
with the author’s elaboration based on the formulae (4) was 
applied. 

However in this case the Nyquist’s frequency limitation is 
also present. The problem is known in the approximation theory 
of discrete sets by using a set of selected functions where the 
functions are required not only to be mutually orthogonal in 
general, but also to be orthogonal on the approximated data set. 
It concerns arbitrary functions , also trigonometric ones. 

The idea of the method of including - to approximation 
model - the function which most contributes to decreasing 
the sum of squares of deviations, is known in a general form, 
and as in the case of invention of the FFT it is connected with 
development and applications of computers [6]. Sometimes 
more firm results may be obtained by using the method of 
rejecting low-significant functions (terms) . In such case the 
approximation is performed in successive steps with the use 
of a model which contains a redundant number of terms. Other 
strategies are also known [6]. In the considered examples of 
run analysis both approaches were applied. For relatively small 
sets of spectral lines (up to 50) a kind of algorithm and proces- 
sing program applied to determining significant spectral lines 
is rather unimportant if confronted with vast computational 
capability of contemporary laptop computers. 

Results of the determination of significant spectral lines in 
the spectrum of the measurement series M with the use of the 
least-squares approximation are presented in Tab.1. 


Tab. 1. Results of the determination of significant spectral lines 
for the run M with the use of TLSA method : q — number of spectral line, 
v —frequency of spectral line, M_ — amplitude, ® — phase, ô — degree 
q q RE. ET 
of decrease of the sum of squares of deviations after including 
the spectral line q into the base of approximating functions . 


In the above presented example the searching for significant 
spectral lines was realized by hand without using the FFT in 
advance. The values v given in Tab.1 were rounded to one digit 
after the point. Some of them were determined with a greater 
accuracy (up to 3 significant digits). Sequence of including the 
spectral lines is of no influence on final result of approximation 
and precise positioning the spectral lines regarding their share 
in decreasing the minimum of sum of squares of deviations is 
not performed in practice but some influence levels are assumed 
instead, which makes the process of searching for spectral lines 
much faster and its automation possible. 

On the basis of the constant and first 10 spectral lines, 
was built the test run M, (Fig.2) which constitutes the least- 
-squares approximation of the run M containing the constant 
and strictly the first 10 amplitude spectral lines (Fig.2) and 10 
respective phases. 


To check the method of searching for significant spectral 
lines the test runs M, , M. and M-p were analyzed also by using 
the TLSA method together with an algorithm for including and 
excluding the spectral lines (Tab.2). 


Tab. 2. Results of the determination of main spectral lines in spectra 
of the test runs with the use of TLSA method : QM,,— spectrum of the run 
M, QM,,, — spectrum of the run M_,,, QM p — spectrum of the run Mp- 


From Tab.2 it results that among the first 10 spectral lines 
for each of the analyzed runs the spectral lines of the same 
frequencies can be found. The obtained result of 100% decrease 
of the sum of squares for 10th spectral line in the spectra QM. 
and QM., is obvious in this case. 

Also, in the case of the run M., (Tab.2) within the set of 
significant spectral lines 10 of them have the same frequen- 
cies, that confirms effectiveness of the method. The level of 
disturbances introduced into the run M_, was significant as the 
maximum amplitudes of the disturbances reached 30% of the 
maximum values of amplitudes of the run M.. 

The absence of correct order of some spectral lines, ob- 
served in Tab.2, as regards the decreasing of the sum of squares 
of deviations is not a mistake but manifestation of features of 
orthogonalization process after inclusion of successive spec- 
tral lines. The analogous result is obtained in every attempt to 
changing this sequence of the spectral lines. Such phenomena 
can happen when values of amplitudes of spectral lines are close 
to each other and their share in approximation of run is similar 
and concerns groups of a few spectral lines (2 or 3) that makes 
determining the spectral lines only a little more difficult. 


SYNCHRONOUS AVERAGING 
AND MOVABLE APPROXIMATION 
OF RUNS 


Run synchronous averaging 


Sampling the torque run M was timed by pulses from shaft 
angle position and for successive rotations of the shaft the an- 
gular positions of sampling points could be assumed identical. 
For each of the shaft rotations 30 measurements initiated by 
shaft angular position, were made. 

Therefore for steady loading conditions of propulsion sys- 
tem the measurements performed within a given shaft rotation 
constitute successive realizations of the same process and the 
deviations from mean values in successive points can be con- 
sidered random ones which have the sum of deviations from the 
mean tending to zero along with increasing number of samples. 
Hence it is reasonable to perform the synchronous averaging of 
run (called also : coherent, linear, prior-detecting or vectorial 
averaging) [2]. It should be mentioned that the averaging is 
a point-wise least-squares approximation procedure. 


Fig.8 shows which way the differences between the mea- 
sured run and its averaged values go depending on a number 
of synchronous averaging cycles. 
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Fig. 8. Runs of the cyclic differences AM = AM; - AM, for some selected 
numbers of synchronous averaging cycles : AM „ — without averaging, 
AM, — for 5-cycle averaging, AM „~ for 10-cycle averaging. 


The runs AM for the whole range of n were achieved by 
extending the run M by the successive run M. To assess effec- 
tiveness of the averaging procedure also the runs of Pearson’s 
correlation coefficient values synchronically calculated can 
be applied (Fig.9). 
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Fig. 9. Runs of the Pearson’ correlation coefficient r,, determined for the 
neighbouring cycles of the run M,,: r „o~ without averaging, t cs ~ JOr 
averaging after 5 cycles, r yc ~ Jor averaging after 10 cycles . 

The runs AM are characterized by the deviations between 
the points distant by one cycle, but the runs r,, are characterized 
by their trends. To determine r,, values for each of them 5 suc- 
cessive samples were taken into consideration. 

In the case of the analyzed data set the maximum possible 
number of averaging was 16. In Fig.10 a fragment of the mea- 
sured run M and the averaged run M, (for 16-cycle averaging) 
is compared, as well as the run of the standard deviations S,, 
is presented. 
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Fig. 10. Comparison of a fragment of the measured run M and the averaged 
run Mo, (for 16-cycle averaging): S,,— standard deviation of the run M . 


C16 
It is easy to observe that the waving period of the standard 
deviations S,, (Fig.10) amounts to about 0.7 of the shaft rotation 
period, hence it can be concluded that this it the large amplitude 
of torsional vibrations of the propulsion system (i.e that fitted 
with one-cylinder combustion engine) which constitutes the 
main cause of torque measurement disturbances. The smooth- 
ness disturbances observed in the run M_,,, may result from 
manufacturing and assembling inaccuracies of the discs. 
The FFT of the single averaged run obtained contains spec- 
tral lines solely for integral values of v (Fig.11). 
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Fig. 11. The amplitude spectrum FM_,, of the run Mo, 


obtained with the use of Matlab software . 
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In Fig.12 a fragment of the averaged run M_..,, was com- 
pared with its two trigonometric approximations for an assu- 
med number of spectral lines, obtained by means of the least 
squares method. 
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Fig. 12. Comparison of the averaged run M, 


and its approximations : 
M,,- that based on the spectral lines q = 1, 2; My- that based 
on the spectral lines q = 1, 2, 3 (Tab.4); Mj,» M3; — their first derivatives . 
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In order to better demonstrate differences between the 
runs Mo, and Mo, the runs of their first derivatives M¿, and 
M}, are also presented (Fig.12). However in both the cases 
shown in Fig.12 one may doubt, making only the use of gene- 
ral knowledge on physical nature of the run M, whether the 
approximation models are correct. Possible adding successive 
spectral lines would decrease quality of the approximation. It 
will be further demonstrated that the problem of smoothing the 
run does not consists in adding or removing the spectral lines 
which appeared in the spectrum of the run M. 

The information removed from the run M as a result of aver- 
aging has been now stored in the run of deviations, AM_,,, whose 
FFT obtained by means of the Matlab, is shown in Fig. 13. 
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Fig. 13. A fragment of the amplitude spectrum FAM „,, obtained 
by means of the Matlab for the run of deviations, AM- = M—M.,,- 


In Tab.3 are presented results of searching for the most 
significant spectral lines in the run AM, , by using the least 
squares method. 
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Tab. 3. Results of determining the most significant spectral lines 
Jor the run of deviations, AM_.,,, by using the TLSA method . 
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In Tab.3 the spectral line of v, = 1 was left in order to 
demonstrate effectiveness of the performed averaging. The 
value of the spectral line amplitude AM, = 0.17 appeared 
significantly smaller as compared with that of the averaged 
run M, = 10.46 (Fig.2, that confirms good effectiveness of the 
performed averaging. 

In addition, the spectrum QAM „1s graphically presented 
(Fig.14) to facilitate making comparisons with the spectrum 
FAM, , (Fig. 13). 
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Fig. 14. Main spectral lines of the amplitude spectrum QAM, obtained 


by using the TLSA method for the run of deviations, AM..,, = M—M.,,- 


The run of deviations, AM..,,, and runs of its 
selected approximations are shown in Fig.15. 
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Fig. 15. The run of deviations, AM_.,, = M - M_,,, and its selected LSA 


approximations: AM,,— that containing the spectral lines for q = 2 - 5, 
AMn- that containing all 10 spectral lines (Fig. 14). 


From the comparison of the runs shown in Fig.15 it results 
that four neighbouring spectral lines of k = 2-5 are mainly 
responsible for the observed ,,beating” (Fig. 14). 

The run of deviations, AM, , can be also subjected to 
further synchronous averaging for selected non-integral fre- 
quencies v and such operation can be repeated in order to split 
the disturbances regarding their characteristic frequencies — in 
an analogous way as it takes place in the case of approxima- 
tion itself. 


Examples of synchronous averaging the run 
of deviations, AM, » for two values 
of v are shown in Fig. 16. 
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Fig. 16. Results of. synchronous averaging of the run of deviations , AM..,, 
of 10 averaging operations, for the various values of v : M z~ a fragment 


of the averaged run for v = 11/30, M „~ a fragment of the averaged 


run for v = 42/30. 


The above mentioned procedures could be applied again 
to the obtained runs or those of differences resulting from the 
averaging in order to finally extract interesting periodic compo- 
nents. Limited accuracy of v - values make a problem in their 
setting. In the case of the run in question they arise from the 
number of samples derived from one shaft revolution. 


The synchronous averaging method should be considered as 
basic one for processing the data of periodical character, mea- 
sured in steady operation conditions of propulsion system. 


Smoothing the runs by means 
of movable approximation method 


In the systems for measuring torque or angular speed the 
main aim of the processing of measurement results is to deter- 
mine in real time the runs of the quantities, characterized by an 
appropriately low level of disturbances. Examples of smoothing 
the run M_.,, with the use of two different approximating objects 
elaborated by this author [9], are shown in Fig.17. 
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Fig. 17. Comparison of the averaged torque run M „and its approxima- 


tions obtained by using the movable approximation method: M „ — the result 
of the triple approximation by using 3rd order power polynomial (Savitzky s 


— Golay 5 filter) within the interval of 11 samples; M_,,— the result 


of the triple approximation by using 3rd order power polynomial 
within the interval of 21 samples and 2 nodes having free constraints 


(cut derivatives: 1" ; 2" and 3"); M}, , M?.,,— 1" derivatives . 


As results from the comparison of the runs M,, and M_.,, 
as well as their derivatives M¿, and M¿,, , the achieved results 
are identical in the applied scale of imaging. The differences 
are negligible and they can be disclosed e.g. by comparing the 
higher order derivatives or remote significant places of values 
of spectral lines of Fourier transforms. As already mentioned 
there is a justified supposition that the approximations by using 
trigonometric polynomials (Fig.12) are loaded by large errors 
resulting from a trigonometric model inadequacy, that can be 
concluded from the comparison of the runs of 1* derivatives 
for both the approximation cases (Fig.12 and 17). The problem 
can be solved by analyzing the motion models of propulsion 
system, e.g. by taking into consideration a run of indicated 
pressure. 

If the Fourier transform for the runs M,, and M_..,, (Fig.17) 
is elaborated it will have the form shown in Fig.18, whereas 
the FFT of the run M, consists 2 spectral lines, and that of 
the run M,, - 3 ones (plus constant) taken from the transform 


FM..,, (Fig. 11). 
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Fig. 18. Ten main spectral lines of the amplitude spectrum QM 


determined by using the least squares method for the smoothed run M, adi 
The transform QM..,, (Fig.18) is characteristic by having 
a dozen or so spectral lines of similarly small values of am- 
plitudes but significantly important for approximation quality 
(smoothing) of the run, which is manifested by a large degree 
of decreasing the sum of squares of deviations after adding 
successive spectral lines (Tab.4). 


Tab. 4. Values of input data to determination process of main spectral lines 
and their parameters for the spectrum QM_.,, (Fig. 18) . 
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In order to fully recognize differences between the runs their 
phase spectra should be also compared. In Fig.19 is shown the 
run of the deviations, AM between the run AM_,, and the 
smoothed run M.,,,. 
: AMeas3 
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Fig. 19. The run of deviations : AM „= Mois- Mea; 


The transform of the above mentioned deviations amounts 
approximately to the difference of the spectrum FM_,,, (Fig.11) 
and QM., (Fig. 18). 

Smoothing by means of the movable approximation method 
can be also applied to the runs of deviations AM, 

Smoothing by means of the movable approximation method 
and decomposition of disturbances can be applied both before 
and after synchronous averaging. A very important advantage 
of the movable approximation is its possible application to 
smoothing and decomposing the disturbances measured in 
unsteady conditions : starting-up, stopping, variable loading. 
The disturbances extracted as a result of decomposition can be 
next subjected to frequency analysis, if necessary. 

Alternative smoothing the runs by applying the Fourier 
inverse transform would require prior multiplication of the 
transforms by appropriate weighting functions. In order to find 
such weighting functions it would be necessary to smooth the 
runs in advance e.g. by applying the movable approximation. 


Interim averaging per cycle 


The interim averaging per cycle is a peculiar case of the 
movable approximation. It is aimed at imaging the variability 
of mean values usually determined as a function of time or 
shaft rotation angle. Such characteristics may be important 
for diagnostics of propulsion systems. In Fig.20 are presented 
results of interim averaging the runs of torque and angular 
speed oscillations, shown in Fig.1. 
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Fig. 20. Results of interim averaging the runs M and œp, (Fig. 1): Mo, ope 
— results of averaging per one cycle (rotation), M — mean value of the set . 


The achieved results can be used in successive steps of 
their analyzing with the use of the above presented methods 
of data processing. 
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CONCLUSIONS 


The methods of synchronous averaging and movable 
approximation should be deemed basic for processing the 
discretely measured runs of torque and angular speed. 


The spectral analysis with the use of least-squares approxi- 
mation of random disturbances and uneven sampling makes 
it possible to extract significant spectral lines in the situation 
of truncated run, random disturbances and non-uniform 
sampling. It also makes it possible to precisely determine 
characteristic frequencies as well as extract significant 
spectral lines of neighbouring frequencies. 


The FFT can be used for preliminary recognizing the spec- 
trum if there is no computer program available for automatic 
determining the significant spectral lines by means of least- 
-squares approximation. 
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a,» b, — trigonometric polynomial coefficients 
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— trigonometric polynomial constant 

— power polynomial coefficients 

— Fourier transform of a given run 

— spectral line number 

— number of trigonometric polynomial coefficients 

— torque, measured torque 

— torque run synchronically averaged per 16 cycles 

— torque test run containing 10 spectral lines 

— test run obtained from the run M, by removing 10 samples 
from each end of set of samples 

— test run obtained from the run M, by adding to it random 
disturbances from the interval <-5, 5> 

— test run obtained from the run M, by its extending by 
zeros as to reach N = 512 samples 

— number of a sample in measurement set : 1, ..., N 

— number of samples in measurement set 

— spectrum of a given run, obtained by using LSA method 

— angular axis, rad 

period 

— number of power polynomial coefficients 

— run of deviations of the run M (a kind of deviations is 
indicated by its indices) 

— run of deviations, AM., = M - Mois 

— degree of decreasing the sum of squares of deviations 
after including the component q to the basis of 
approximating functions 

— frequency 

— dimensionless spectral line frequency, v, > 0 

— time 

— phase. 


Indices 


— number of a distinguished (significant) component 
of a spectrum (spectral line) in Fourier transform 
— limit value 
— successive number of spectral line in Fourier transform 
— number of spectral line in the transform obtained 
by means of TLSA method 
— symbol of 1* derivative. 


Abbreviations 


Discrete Fourier Transform 

Fast Fourier Transform 

approximation by means of trigonometric polynomials 
and least squares method. 
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on Exploitation 
Foundations 


On 23 March 2006 Faculty of Ocean Engineering and 
Ship Technology, Gdansk University of Technology, ho- 
sted the members of the Regional Group of the Section on 
Exploitation Foundations, Machine Building Committee, 
Polish Academy of Sciences ( PAS), who took part in their 
successive scientific seminar. 


The seminar’s scientific program consisted of four pa- 
pers presented by the scientific workers of the Faculty : 


= Comparative analysis of main engine loads and output 
of multi-bucket dredgers by using probabilistic appro- 
ach — by D. Bocheński 

= Diagnostics of the energy transformation process 
of steam turbine power plant of large output — by 
J. Głuch 

= Research on water-lubricated short-sleeve bearings of 
propeller shaft — by W. Litwin 

= A concept of application of artificial neural networks to 
represent reliability probabilistic models — by J. Rud- 
nicki 
After interesting discussion on the presented topics the 

seminar participants were acquainted with the Faculty’s 


Pi 
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Calculation of the mean long-term 
service speed of transport ship 


Part I 
Resistance of ship sailing on regular shipping route 
in real weather conditions 


Tadeusz Szelangiewicz 
Katarzyna Zelazny 
Szczecin University of Technology 


ABSTRACT 


Service speed obtainable by a ship in real weather conditions when sailing on a given shipping route, is 
one of the major parameters which have great impact on ship operation costs. The so far used, very ap- 


ï 


proximate method of service speed prediction based on “service margin”, is very little exact. In this paper 

a new method based on additional ship resistance dependent on mean statistical parameters of wave and 

wind occurring on a given shipping route, is presented. The mean long-term service speed is calculated on 

the basis of the calculated additional resistance and the screw propeller and propulsion engine parameters. 

Also, a new definition of service margin and a way of its calculation is presented apart from the results of 
the mean service speed calculation depending on ship 8 type and size and shipping route. 


Keywords : ship service speed, wind, waving, shipping route, service margin, long-term prediction. 


INTRODUCTION 


The ship owner, when planning building a new ship, on 
the basis of his experience and knowledge of shipping market, 
determines the ship’s deadweight and predicts its service speed 
to be reached by the ship on a given shipping route. The service 
speed obtainable by a ship in real weather conditions (mainly 
wind and waves) is one of the most important parameters 
greatly influencing ship’s profitability on a given shipping 
route. In the design algorithms used today in ship design offices, 
ship’s deadweight is determined with a satisfactory exactness, 
but calculation of ship’s service speed is rather inexact and 
practically known already after ship’s delivery and some time 
of its service. Hence in the ship building contract, is given 
another speed, so called contract speed which the ship has to 
achieve in still-water conditions. On the basis of still-water 
resistance characteristics and rated output of propulsion engine, 
the obtainable ship service speed is predicted by accounting 
for relevant coefficients, Fig.1. 

In that prediction an increase of ship resistance in service 
conditions in relation to that in still water conditions (rep- 
resented by the service margin SM in Fig.1) is not exactly 
determined as it depends on real weather conditions occurring 
on a given shipping route. 

Basing on known statistical data on winds and waves 
possible to occur on a given shipping route, one is capable in 
determining a long-term distribution function of additional 
resistance, and next, for an assumed main engine output, 
a long-term distribution function of ship’s speed possible to 


be obtained on a given shipping route. Such approach makes 
it possible not only to exactly determine the mean statistical 
service speed of the ship, but also to investigate the influence 
of power output of a selected main engine on the ship’ s service 
parameters including its long-term service costs. 


A 


Propulsion power 


VE V 


Ship speed 


Fig. 1. Predicted service speed of a ship. Notation : Nn — rated output of 
engine, PS — shaft power, PD — power at propeller s cone, OM — operatio- 
nal margin, SM — sea margin, V,,— contract speed, V,, — predicted service 
speed, 1 — propeller curve for still-water conditions, for non-fouled hull, 
2 — predicted propeller curve with service margin, for real conditions . 
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Attempts to determine the service speed of ship in real 
weather conditions have been already reported in [4, 23, 24], 
however there were not accounted for statistical distributions 
of wind and wave parameters met on a given shipping route. 


STATISTICAL LONG-TERM 
PARAMETERS OF WINDS AND WAVES 
ON SHIPPING ROUTES 


In order to determine the mean long-term ship service 
speed the ship’s resistance in real weather conditions should 
be determined in advance. The statistical long-term wave 
parameters, namely : H, - the significant wave height, T, - the 
mean characteristic period, u - the mean geographical direc- 


tion, were determined on the basis of multi-year measurements. 
Results of the measurements are contained in the atlas of wave 
parameters [10] and in extended form in [11], in which seas 
and oceans have been divided additionally into the set of sea 
areas shown in Fig.2 ( where each sea area is characterized by 
its own statistical wave parameters). 

The atlases [10] and [11] contain sample size of waves char- 
acterized by H, and T, parameters and recorded in a given area 
for a given season and wave direction u, as exemplified in Tab. 1. 
On the basis of these data, was elaborated a computer database on 
the waves met on seas and oceans, which contains the occurrence 
probabilities of waves with H,, T, parameters, f,,, , and the oc- 
currence probabilities of u directions , fa , for all the sea areas 
(split into seasonal data), as exemplified in Tab.2 and 3. 


170 60 i50 MO BO RO IIO IOO 90 80 70 60 5O 40 30 20 0 O 1 29-30 40 50 6O 70 80 90 100 lIỌ 20 130 MO ISO 160 17% i80 170 


70 1 RRR Seal NENA ES EE as E 
170 6O ISO 140 130 [20 IIO 100 90 BO 70 60 SO 40 30 20 10 O 1 20 30 40 50 60 70 80 90 IOO IIO 120 130 MO I5Ọ 160 170 180 170 


70 


Fig. 2. Oceans divided into set of sea areas, acc. [10] . 


Tab. 1. Sample size of waves of H, height and T, period for a given w direction and sea area [10] . 


Sea area: 1 


Season of the year : December — February (winter) 


Wave direction : u = 000 deg 


10+ 11 


12+13 | 14+15 
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Tab. 2. The occurrence probabilities , f,,, , of wave of H, height and T, period for a given \ direction and sea area [10] . 


Sea area: 1 


ONOMY 


Y 


Season of the year : December — February (winter) 


Wave direction : u = 000 deg 


0.091 | 0.023 | 0.011 | 0.011 


| | 
0.011 i 


0.011 | 0.045 | 0.057 


0.023 


0.023 | 0.011 


0.011 | 0.011 | 0.045 | 0.011 | 0.011 


-— [oon | 
Ee ae 
| [oon | 


Tab. 3. The combined probability, f,,,° f, , of occurrence of wave of H, height and T, period, 
with taking into account the probability f£ for each direction u in a given sea area, [10] . 


Sea area: 1 


Season of the year : December — February (winter) 


Wave direction : all 


0.0009 


0.0026 | 0.0094 | 0.0099 | 0.0034 | 0.0030 | 0.0013 


0.0094 | 0.0081 | 0.0069 | 0.0039 | | 
0.0009 | 0.0051 | 0.0017 | 0.0009 | | 0.0004 


0.0004 


O 
| | 
| 
[0.0009] | 
| 
| 
a 


| 75 | 0.0017 | | 0.0004 | 0.0017 | 0.0073 | 0.0064 | 0.0013 | 0.0004 | | [| 0.0004 
| 80 [0.0013 | 0.0004 | 0.0009 | 0.0009 | 0.0056 | 0.0030 | 0.0004 | 0.0009 [ [| | 
| 85 | | [0.0004 | 0.0013 [0.0038 | 0.0017 | 0.0004 | 0.0004] Ù Ù |] 


| [| 0.0004 | 0.0013 | 0.0034 | 0.0004 | 0.0004 eS as 
| [0.0026 | 0.0026 | 0.0034 | 0.0009 | 0.0013 | 0.0021 | | 0.0004 


It 
0.0009 | | 


| [0.0004 


pS 0 PT 


The statistical wave parameters used for calculation of ship 
resistance due to wind and waves should be correlated with the 
wind parameters: V, - mean wind speed, y, - mean statistical 
geographical direction of wind. As such atlases for wind do not 
exist the appropriate relations between wind and wave parame- 
ters were elaborated on the basis of measurements [3 , 25]. 


Basing on literature sources, one elaborated the relations V, = 
= V (H, , T,) for the purposes of this work, assuming that : 


+ Wind direction coincides with wave direction : 
* 


Ya HU (1) 


* 


+ The occurrence probability of wind direction, f,, , 
is the same as that of wave direction , fi i 


f 


yA fı 


(2) 


+ The occurrence probability of mean wind speed, f,, , 


f . 


is the same as that of wave of H, , T, parameters , f,,, : 


fin E For 
Results of calculations of the relations V, = V,(H, , T,) 
for the wave parameters contained in [10], 
were published in [25]. 


in the literature on wind-generated sea waves it is assumed that if the difference between wind and wave directions does not exceed 30° 


then both the directions are the same; such state occurs during 80% time of observations on seas and oceans [6, 22] 
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SHIP RESISTANCE 
IN REAL WEATHER CONDITIONS 


In real weather conditions during sea voyage of a ship on 
a given shipping route, its total resistance consists of still-water 
resistance and additional components : 


Re =R+AR (4) 
where : 


R -—still-water resistance of ship (accounting for the drift 
angle and possible surface sea current) 

AR — additional ship resistance resulting from real sailing 
conditions : 


AR=K 7 KK. (5) 


R., — additional wind-generated resistance 

R — additional wave-generated resistance 

R,, — additional resistance resulting from e.g. rudder 
action to keep a given ship course. 


STILL-WATER RESISTANCE OF SHIP 


The ship’s still-water resistance is usually measured (by 
model testing) or calculated for the ship in rectilinear motion. 
In real weather conditions the ship sails under a certain drift 
angle due to oblique action of wind and waves or/and possible 
surface sea current. Hence the still-water response for the ship 
at steady speed, is composed of : 


1 


R,= zP WC 
1 2 

Ry= FP SVavCy(Brv) (6) 
1 

M,= 3P SiViy C..(Bry) 


where : 


R,, R, , M, — components of still-water resistance forces and 
moment of ship sailing at drift angle B,.,, and 
surface sea current, respectively, (Fig.3) 

p  — water density 
S — lateral projection of underwater ship hull 
surface onto ship’s plane of symmetry (PS) 
Vœ — relative ship speed 
Pey ~ relative drift angle 
L  — ship length 
Cis C, , C,, — coefficients of resistance forces and moment. 


xo ft 


In Eqs. (6) the relative speed 
and relative drift angle is as follows : 


Vrv = View + a (7) 
-V 
Bry = arctg — (8) 
RVx 


Vavx = We ~ Ve cos Be 
| (9) 


where : 
V -— absolute ship speed 
V = Vcos B \_ absolute ship speed 
V,=—Vsin B components 
B — drift angle at absolute speed 
V- surface sea current velocity 
Pe — sea current direction relative to ship 


Be= Ve T Y (10) 

Ye —~ geographical direction of surface sea current, (y, = 

= 0° northbound current, y, = 90° eastbound current) 

y — geographical course of ship, (y = 0° northward course, 
y = 90° eastward course). 


If the sea current velocity V, = 0 then the absolute ship 
speed V and absolute drift angle B is valid for Eqs. (6). 


In the subject-matter literature can be found various em- 
pirical formulae based on the tests of many ship models, by 
means of which values of the coefficients C., C, Ch (Bry) 
can be calculated. In [20], the coefficients C, , C, , Ca for the 
directions defined in Fig.3, have the following form : 


Xo 


075 
C, =0.5C, sin 2Bpy COSBpgy + 


j- Prv 


k 


X 


C, =0.075sin 


180° — arcsin 
0 


+C,sin*Bpy+C3sin* 2Bpy (11) 


Cn =m, sin 2Bgy +m, sin Bgy + 


+m; sin? 2Bpy+m,sin* 2Bpy 
where : 

Cas k C, C, C, , m; , m, , m, , m, — the coefficients 
dependent on ship hull parameters, and S = L- T : o (where: 
T — ship draught, o — reduction factor for underwater hull 
surface side area) is the surface used in this method. 


00 


Yo 


Fig. 3. Relative ship speed and drift angle in still water with surface sea current . 
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An example of calculated values of the coefficients C, , 
C Cn (Bey) acc. [7] is shown in Fig.4. The shown coefficients 
were determined from the model tests [18] but only for small 
values of the drift angle By (where, as often used, the angles 
are linearly dependent on the angle By ). 
0.04 4 
0.03 4 
0.02 4 
_ 0.01 | 
= 0 


O -0.01 4 20 40 60 80 


pre] 


T T T T id 


100. 120 140 160 180 


-0.02 4 
-0.03 4 


0 20 40 60 80 100 120 140 160 180 
BI] 
v 20 40 60 80 100.120 140 160180 
-0.05 5 ge 
-0.14 


----- model tests [18] calculations [19] 


Fig. 4. Ship resistance coefficients for the ship K1* during motion at the 
speed V = 8.4 m/s and a drift angle, as calculated, acc. [20] . 


Inserting the relations between the drift angle By 
and the ship speed V,,, into Eq. (11) one obtains : 


arccos Vv 
C, = 0.075sin||180°— arcsin—*_] 1- RY 
0.07 k, 
7 VivsVevy Vaevy|Vevs 
C==G-—__, --&- 3 +t 
Vev Vev 
VV, (12) 
-16C RV. 2 a 
RV 
C,,=—2m, NevsVivy —m, Vavy + 
RV Vav 
y? y? V: NV: V. 
Z 8m, — RVy ~16m, RV: Rw ie 
RV RV 


* the parameters and dimensions of the ships for which the relevant calcula- 


tions will be presented in Part III of this paper. 


The coefficients C, C, C, in the form (12) make it possible 
to easily calculate the drifting speed Vi and thus the drift angle, 
provided the components of the external forces exciting ship 
drift are known. 

The coefficient Cy which appears in the first of Eqs (12), can 
be calculated from still-water resistance for ship in rectilinear 
motion and without drift : 


c= Ca) 


Xo 1 
ar> 


where : 


(13) 


C,,(V) — the still-water resistance coefficient 
for ship in rectilinear motion : 


R 
CR(V) = v2 


and 
R - still-water resistance for ship in rectilinear motion. 


(14) 


As the coefficient C,(V) has no constant value (Fig. 5) hence 
the resistance R or coefficient C, will be provided for a given 
ship also in the form of the table of discrete values dependent 
on the speed V. For a given ship its resistance can be approxi- 
mately measured by means of model testing or derived from 
approximate formulae e.g. those provided by Holtrop-Mennen 
method [13, 14, 15] or Hollenbach method [12]. 

The resistance R to a large extent depends on hull surface 
state. The above mentioned methods make it possible to take 
into account a real state of hull surface and thus to investigate 
its impact on ship’s service speed, that has been used in the 
presented method. 

X 


500 + R [KN] 


v 


400 + 


300 + 


200 5 


100 + 


V [m/s] 


1 3 5 7 9 


Fig. 5. Example characteristics R(V) and C,(V) for the ship K1 . 
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WIND INFLUENCE ON SHIP IN MOTION 


The mean wind forces acting on the ship in motion 
can be calculated by using the following formulae : 


1 
Rya= —5PaSx VeaCax(Bra) 


1 
Rya= z PaSy Vaal Pra) (15) 


1 


M,a = 5 PASy LV RAC AnCBra) 


where : 


Pa — air density 

S,» S, — areas of front and side projections of above water 
part of ship onto midship and symmetry plane of 
ship, respectively 

— ship length 
— relative wind speed (Fig.6) 

Cy > Cam (Pra) ~ aero-dynamical drag coefficients of 
the above- water part of ship surface, dependent 
on the relative wind direction (B, ,) 

— relative wind direction (Fig.6) 


2 2 
Vka = 4 Vax t VRay 


L 
Vaa 
Ce 


Bra 
(16) 


l (17) 
Vray= Va Sin Ba 
B=), Ty + 180° (18) 


where : B, — wind direction relative to ship 


-V 
Ry (19) 
RAx 


V — absolute wind speed 

Ya — geographical direction of wind (y, = 0° 
— north wind, y, = 90° — east wind) 

V — ship speed 

y — geographical course of ship. 


Bra = arctg 


In Eqs. (15) the ship drift angle B was omitted because of its 
small value and very low influence on the relative wind speed 


Bas and thus on the wind impact forces. 
x 
of 


0o 


Fig. 6. The reference frames of coordinates, 
speeds and directions of ship and wind . 


28 POLISH MARITIME RESEARCH, No 4/2006 


One can calculate the wind forces acting on the ship in 
motion at a given ship speed V and course y, assuming the 
wind parameters (V, , y,), the aero-dynamical drag coeffi- 
cients Ci,» Ca > Ca calculated by means of the approximate 
formulae. given in like [21] or derived from the model tests 
[2] for the ship K1 as shown in Fig.7. In Eqs. (15) the force 
R, represents the additional resistance due to wind, and the 
force R,, and moment M,, may induce ship drift and change 
of its course, and ship’s steering gear must be used to balance 
it, and due to the gear action the additional resistance R,, will 
be then produced. 


0.8 + 
0.6 4 
0.4 4 
0.2 4 


T T T iai 


120 140 160 180 


PraL’] 
0 ba T T T T T T T T T > 
100 


120 140 160 180 


100 120 140 160 180 


Fig. 7. The aero-dynamical drag coefficients C, Cy Ca. 


for the ship KA, acc. [2] . 
WAVE INFLUENCE ON SHIP IN MOTION 


The mean irregular-wave—generated forces (i.e. mean wa- 
ve-generated drift forces) acting onto ship in motion can be 
calculated by using the formulae : 


Be 
Ry = 2pe——|Cwx(/By.V)S_c(@)de 
0 


Bee 
Ryw = 2pg——|Cuy(o/By.V)S:c(o)do (20) 
0 


Mw = 2pgB7)Cwn(@/BwsV)Sec(@)do 
10) 


where : 


— water density 
— acceleration of gravity 
— ship breadth 
v Cw , Cym (@ / By» V) — coefficients of regular- wave- 
-generated drift force dependent on the wave di- 
rection relative to ship, B,, and ship speed V 


Q Wn Dd 


0) — regular wave frequency 
B,, ~ wave direction relative to ship (Fig.8) 
— = fe} 
B, =H- w+ 180 (21) 
u — geographical direction of waves, (u = 0° - north 


wave, u = 90° - east wave) 
S.(@) — wave energy spectral density function (dependent 
i on the significant wave height H, and mean wave 
period T). 


e 


Oo Yo 


Fig. 8. Mean wave forces acting onto ship . 


Assuming the wave parameters (H, , T,, 11), and the coef- 
ficients of wave — generated drifting force, C,, Cw > Cwm > 
(calculated acc. [1 , 5] or derived from the model tests as for 
the ship K1, shown in Fig.9), one can calculate, for a given 
ship speed V and course w, the wave-generated forces acting 
on the ship in motion. As in the case of wind action the force 
R w represents the additional ship resistance due to wave ac- 
tion, and the force Rw and moment M „ may cause ship drift 
and change its course. 


ADDITIONAL SHIP RESISTANCE 
DUE TO PASSIVE RUDDER 


When ship sails in waves, especially when oblique wind 
and waves influence the ship’s motion, are generated lateral 
forces and moments which force the ship’s course changing 


A 
Xo 


V = 0 [m/s] 


By = 180° 


Fig. 9. Wave-generated drift force coefficients C 


for the ship K1, acc. [8]: 


CoG 


Wx ? “Wy ? “Wm 


a) for various wave directions relative to ship at V = const, 
b) for various ship speed at B, = const . 


Fig. 10. Forces acting on passive rudder . 
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and drifting. In order to keep the ship’s course constant the 
rudder blade must be inclined, that produces the additional 
resistance R p- 

In the literature on ship manoeuvring several algorithms 
for calculation of hydro-dynamical forces on passive rudder, 
including those dealing with additional resistance , can be 
found, e.g. in [9 , 16]. According to [16], the passive rudder 
forces can be calculated by using the formulae : 


RR = 


M g = a,Fy cos Op 


Fy sin 8,| 


ayFycosdp (22) 


where : 


6, -passive rudder angle 
a — coefficient of hull influence on the rudder force R. 


a = coefficient of hull influence on the rudder moment 
Mrz 
a= aX, (23) 


X, — abscissa of rudder axis measured from the ship mass 
centre G (X, < 0) 
F  — rudder normal force 


1 613A 
NOP A4225 


à — rudder aspect ratio 
A, ~ rudder surface area 
Vp ~ water inflow velocity to rudder 


a, — effective rudder angle of attack. 


AkVŽsin apg (24) 


As a result of passive rudder inclination the moment Mp 
of force R „ appears, and in order to keep the set course of the 
ship, the rudder moment should have such value as to balance 
the resultant forcing moment due to action of wind, water flow 
and waves: 


M, =— M, (25) 
where : 
M,r = M,a + M,w + M, (26) 


Hence, value of the rudder angle 6,, calculated by using 
Eq. (22) under assumption of keeping the ship course constant, 
will be dependent on the wind and wave parameters and ship’s 
drift angle. 


NOMENCLATURE 


A — rudder surface area 
— coefficients of hull influence on rudder forces 
— ship breadth 
— aero-dynamical drag coefficients 
— still-water ship resistance coefficient 
— water resistance coefficients of forces 
and moments 
— coefficients of regular -wave -generated drift 
force 
approximation coefficients dependent on form 
of underwater part of ship hull 
F — rudder normal force 


CoC 


Wx’ ~Wy ? “Wm 


C C, C, C97 


f — occurrence probability of wave of H, and T, 
parameters, coming from u direction 

f — occurrence probability of mean wind speed 

ta — occurrence probability of wind direction 

f, — occurrence probability of wave direction 

g — acceleration of gravity 
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H, — significant wave height 

k, — approximation coefficient dependent on form of 
l underwater part of ship hull 

L — ship length 

Mr — resultant forcing moment due to wind, water 


flow and waves 

— approximation coefficients dependent on form 
of underwater part of ship hull 

R — still-water resistance of ship 

— total resistance of ship 

— still-water ship resistance components of force 
and moment during sailing at drift angle B, 
respectively 

— mean wind forces and moment acting on ship in 
motion (where R, - additional wave-generated 
ship resistance), respectively 

— passive rudder forces and moment (where 
R p- additional rudder-generated resistance), 
respectively 

— mean wave-generated drift forces and moment 
(where R œ- additional wave-generated 
resistance), respectively 

S — area of side projection of underwater part of 

ship hull onto its plane of symmetry 

S — areas of front and side projections of above 
water part of ship onto midship and symmetry 
plane of ship, respectively 

— wave energy spectral density function 

— ship draught 

— mean characteristic wave period 

— ship speed 

— wind speed 

— sea current speed 

— water inflow velocity to rudder 

— relative wind speed 

— relative ship speed 

— ship speed components 

— rudder axis abscissa measured from ship mass 
centre 


FR 


<< << <4 <n 
< 


eee 
< 


— effective rudder angle of attack 

— ship drift angle 

Ba — wind direction relative to ship 

— current direction relative to ship 

— relative wind direction 

— relative drift angle 

— wave direction relative to ship 

Vin — geographical direction of wind 

Yc — geographical direction of sea current 

— additional resistance due to rough weather 
Ù — passive rudder angle 

A — rudder aspect ratio 

u — geographical direction of wave 

p — water density 

Pa — air density 

y — geographical angle of ship course 
(0) — regular wave frequency 
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On 6-7 October 2005 Faculty of Naval Architecture 
and Ocean Engineering, Istambul Technical University 
(ITU) organized : 


8" International Ship Stability Workshop 
which was initiated by the keynote address on : 


Effects of initial bias on the roll response 
and stability of ships in beam seas 
by A.Y. Odabasi and E.Ucer (ITU). 


Program of the scientific meeting consisted 
of 6 topical sessions : 


Theoretical prediction of intact stability (4 papers) 
Experimental investigation of intact stability 

(4 papers) 

Special problems on ship stability (4 papers) 
Theoretical development in damage stability 

(3 papers) 

Assessment of ship stability safety (5 papers) 

as well as the special session on : 


(— > +e ef 
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SHIP STABILITY WORKSHOP 


A Regional issues and activities with one paper on : Anato- 
my of a capsize : then and now — by M. Tyalan (ITU). 


The papers presented during the topical sessions were 
prepared by representatives of the universities and scientific 
research centers from Japan and USA (4 papers each), Greece 
and UK (3 papers each), Poland and Russia (2 papers each), 
Brazil, Holland, Italy, Korea, Norway, Sweden and Turkey 
(1 paper each). 


Polish scientific workers contributed 
in preparation of the following papers : 


© Appraisal of risk assessment approach to stability of ships 
— by L. Kobylinski (Foundation for Safety of Navigation 
and Marine Environment Protection, Poland) 

€ Risk characterization of the requires index R in the 
New Probabilistic Rules for Damage Stability — by M. 
Pawtowski (Academic Visitor to Ship Stability Rese- 
arch Centre (SSRC); Gdansk University of Technology, 
Poland), and D. Vassalos (SSRC; The Universities of 
Glasgow and Strathclyde, UK). e 
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OPERATION & ECONOMY 


This is the name of successive project conducted 

by Faculty of Ocean Engineering and Ship Technology, > 
Gdansk University of Technology 

within the frame of EUREKA European projects 


EUREKA 


PROJECT DESCRIPTION 


Program aim is to introduce into ship structural design a new technology - laser welding techniques. Such techniques make 
possible to manufacture steel sandwich panels - units which are able to replace some ship hull structural parts like bulkheads, 
inner bottom structures, flat side surfaces or internal decks, walls etc. It is expected that application of steel sandwich panels will 
reduce manufacturing costs and also increase internal space. Additionally such ‘twinning’ of shell plating increase resistance of 
structure on perforation which rise safety factor in relation to environmentally hazardous cargoes. Although the basic strength 
characteristics of steel sandwich panels are known as result of previous research programs, implementation such new structure 
into classical ship hull design requires new knowledge. Project will be focused on complex analysis of simplified structure of 
typical small and medium size Baltic operated ship in order to replace some parts by sandwich panels. Detail design of new 
structure especially design of joints will be performed. Numerical calculations and real scale laboratory tests of joints will be 
carried out to verify applied design approach. Analysis of technical and economical aspects of application of laser welding 
techniques in ship manufacturing process will be carried out. 

The expected result is to create a technically advanced and cost effective new generation of small and medium size Baltic 
operated ship with implementation of the economically profitable and environmental friendly new technology. Expected re- 
sults of the project will lead to application of new manufacturing techniques, saving of material and manufacturing costs and 
reduction of weight of ship structure. 


GENERAL AIM on 


Replacement of shell plating with two-directional system E 
of stiffeners by two shells stiffened by internal unidirectional 


laser welded stiffeners system. | F200 x 90 x 9 


Expected results of application 
of the SANDWICH structure : 


= weight reduction: up to 34 % 
= costs reduction: about 50 % 
= increasing of internal space 


ACHIEVEMENTS 


Design of a technically advanced and cost effective new 
generation of small and medium size Baltic operated ship 
structure with implementation of the economically profitable 
and environmental friendly new technology. 

Expected results of the project will lead to application of 
new manufacturing techniques, saving of material and manu- 
facturing costs and reduction of weight of ship structure. Effect 
of such approach will give reduction of environment pollutions 
during manufacturing processes, reduction of manufacturing 
costs due to new technology implemented as well as reduction 
of ship operating costs due to better economic indicators. 

Area of potential application will be mainly shipbuilding 
industry but there are some ideas to spread such approach 
into other market sectors like automotive, railway or civil 
engineering. 
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